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ABSTRACT 


The vertical discretization in a linearized baroclinic prediction model was analyzed 
by comparing various finite element and finite difference solutions following Jordan 
(1985) and Shapiro (1987). The baroclinic instability experiments of Shapiro (1987) were 
augmented to include the unstaggered vertical scheme from Jordan (1985). Two basic 
wind profiles were used and the experiments were run with different resolution models 
and horizontal wavelengths. For a given wind profile and vertical resolution, different 
models performed better. The finite element models for the staggered vertical grids did 
not perform up to their possibilities due to the boundary elements. However, for the 
unstaggered vertical grid, the finite element model did better than the finite difference 


model in most cases. 
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I. INTRODUCTION 


At the present time, most numerical weather prediction models use finite differences 
to accomplish the vertical discretization even though they use finite difference, finite el- 
ement, or spectral horizontal discretizations. The Canadian regional and hemispheric 
models (Staniforth and Daley, 1977 and 1979) which use finite elements in the vertical 
are an exception. The successful numerical prediction of synoptic evolutions requires a 
proper representation of the vertical variation of the predictive fields. Since smaller scale 
features such as fronts (Hoskins and Bretherton, 1972 and Willams, 1967) and the large 
scale planetary waves (Gall, 1976) are forced by energetic svnoptic-scale features, it fol- 
lows that all predictive scales of motion mav be sensitive to the vertical discretization 
used in the numerical mdels. 

Most of the finite difference vertical discretizations use a staggered arrangement of 
variables. Winninghoff (1968), Arakawa and Lamb (1977) and Schoenstadt (1980) 
demonstrated that staggering of variables in the horizontal improves geostrophic ad- 
justment and the response to small scale forcing. Most quasi-geostrophic models 
(Charney and Phillips, 1953) use vertical staggering where the vertical motion and the 
temperature are carried between the levels which carrv horizontal velocity and pressure. 
This arrangement will be referred to as grid B. Lorenz (1960) introduced a different grid 
for the balance equations which was designed to conserve energy. This arrangement 
places only the vertical velocity between the levels which carry the other variables (hor- 
zontal velocity, pressure and temperature) and will be referred to as grid A. Tokioka 
(1978) analyzed a number of vertical grids with linearized equations and found that grid 
A has a computational mode in the temperature field. Arakawa (1984) compared 
baroclinic instability for grids A and B in the linearized quasi-geostrophic equations. 
He found a false short wave instability for grid A which did not occur with grid B. This 
problem is related to the computational mode in the temperature field. Another diffi- 
culty with grid B is that the matrix which must be inverted to find the temperature from 
the pressure 1s singular. This is especially important for initialization. Many operational 
primitive equation models use grid A for energy conservation (Shapiro, 1987). 

The use of finite elements for the vertical discretization can be expected to give a 
more accurate representation of vertical variations. The finite element method ts a spe- 


cial case of the Galerkin procedure which represents the dependent variables with a 


Weighted sum of basis functions that have a prescribed spatial structure. The finite ele- 
ment method emplovs basis functions which are zero except in a limited region where 
they are low-order polynomials. This method has been used in engineering statics (e.g., 
Zienkiewicz, 1977) and it has been applied to fluid dynamics and hydrology (see Gray 
and Pinder, 1976). The finite element method has been successfully applied to 
meteorological prediction with the shallow water equations by Cullen (1973), Hinsman 
(1975) and Staniforth and Mitchell (1977, 1978). Cullen (1973), Neta er al (1986), and 
Neta and Williams (1986) demonstrated that finite element formulations with piecewise 
linear basis functions are more accurate than second order finite differences. 

Jordan (1985) compared six linear, baroclinic, vorticity-divergence equation models 
using three grid schemes, grid A, grid B and an unstaggered grid. The three grids are 
depicted in Fig. 1. Shapiro (1987) improved the models for grids A and B, and corrected 
several problems in the heating and mountain forcing experiments. He also ran a 
baroclinic instability experiment without heating or mountain forcing using linear verti- 
cal shear in the mean wind profile. With the ever-increasing density of observations and 
the advent of satellite-based cloud pictures, transient disturbances of subsynoptic scales 
have been observed in certain preferred locations such as polar lows (Satyamurty ez al, 
1982). Satyamurty er a/ (1982) found that unstable modes can be generated at smaller 
horizontal wavelengths using hyperbolic tangent profiles. Farrell (1985) suggests that 
for backward tlting initial states the waves can experience rapid growth in the early 
stages and then settle into the most unstable growth. 

The purpose of this study is to augment the baroclinic instability study by Shapiro 
(1987) by including the unstaggered models from Jordan (1985) and by looking at 
smaller resolution models and smaller wavelengths. Various hyperbolic tangent vertical 
profiles of the mean wind field are used to focus on more difficult baroclinic instability 
situations with smaller horizontal wavelengths and vertical scales, and with greater sen- 
sitivity to large local shears. This study also experiments with initial states which have 
a backward tilt with height in an effort to investigate the enhanced initial growth due to 


the tilt. The results of the experiments are given in Chapter III. 
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Fig. 1. Three vertical grids. 


Il. NMIODEL DESCRIPTIONS 


A. MODEL FEATURES 

Jordan (1985) developed six numerical models with several features to make easy 
modifications for a wide range of experiments. Shapiro (1987) modified the A and B 
models to improve the performance of the models near the boundary and used numerical 
integration to evaluate the entries of the matrices. He also added menus to make the 
transition between four experiments simple. The user is able to prescribe heating, 
mountain topography, velocity perturbation, or baroclinic experiments and the model 
will make the prescribed changes in the variables governing these cases. The vertical 
velocity vector in finite element model B has been corrected to include an additional el- 
ement in an attempt to eliminate the oscillation in the temperature profile from Shapiro 
(1987). The two models for grid C have been modified to run baroclinic instability ex- 
periments and to include the same menus as the A and B models. Another menu has 
been added to all six models to select between linear and nonlinear wind shear profiles. 
The models are written in modular structure using FORTRAN 77. There 1s parallel 
construction between models. The subroutines used in one model are verv similar to 
those used in the other models. The models can run quickly on an IBM-3033 
mainframe. For example, a 96-hour forecast for a 6-laver finite element model uses less 
than six seconds of computer processing time; however, the same forecast for 60 lavers 


can take up to 220 seconds. 


B. GOVERNING EQUATIONS 

Each model approximates the same set of governing equations. The vorticity 
equation (2.1), the divergence equation (2.2), the surface geopotential equation (2.3), 
and the first law of thermodynamics (2.4) are the prognostic equations for the forecast 
variables vorticity, divergence, surface geopotential and potential temperature. The 
surface geopotential equation is the lower boundary condition on the vertical velocity. 
The vertical coordinate Z = — In(p/p,) is used, but the non-boussinesq terms involving 


e-* are replaced by one. The prognostic equations in the coordinates x, y, Z and t are 


dg 6w dv éw ou 
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In these equations, 
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@ is the geopotential, ¢ = gZ, 

@, is the surface geopotential. 

T is the potential temperature, 

u is the x-component of velocity, 

Vv is the y-component of velocity, 

w is the vertical velocity, 

Q is the diabatic heating per unit time per unit mass, 

MTS is the forced vertical velocity due to flow over mountain topography, 


fis the Coriolis parameter, 


B is df/dy, 
di)  o&) c{ ) c( ) et ) 
fe, | CU ax ap i + W az and 


V? is the horizontal Laplacian operator. 


The prognostic equations are linearized by expanding the variables into their mean 
and perturbation states, as in Jordan (1985). The resulting linearized forecast equations 
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where R is the gas constant for air, (’) denotes perturbation quantities and(__) denotes 
mean quantities. The use of Xbar in the text will be used to denote mean quantities of 
a variable X (Shapiro, 1987). 

The diagnostic variables, u’, v’, w’ and ¢’, are calculated from the forecast variables 
using the definitions of divergence, vorticity, the hydrostatic equation and the continuity 


equation. The relationships are given in equations (2.9) through (2.12). 
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The use of primes to denote perturbation quantities will be discontinued. All quantities 
used in the remainder of the paper will be perturbation quantities unless otherwise noted. 

The mean state is assumed to be in hydrostatic and geostrophic balance. The term 
OTbar/dy in the first law of thermodynamics can be evaluated by taking ¢0/éy of the 
hydrostatic equation and substituting for déd@bar/déy from the geostrophic relation, 
Opbar/éy = —fubar. Thus, 
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Geostrophic balance of the mean state at the surface implies 
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The expressions (2.13) and (2.14) are substituted into equations (2.8) and (2.7), respec- 
tively. 
A singlewave spectral representation 1s used in the x-direction, with wave number 


u=27/L, where L 1s the wavelength in the x-direction. The perturbation quantities have 


the form 
C(x,Z,t) = A,(Z,t) cos ux + A,(Z,t) sin px, 15) 
D(x,Z,t) = D,(Z,t) cos ux + D,(Z,t) sin px, (2.56) 
MixGZ ,t) = 1,(Z,t) cos ux + 1,{Z,t) sin ux, (2.17) 
ket) — S,(Z,t) cos nx + S,(Z,t) sin ux, (als) 
BixeZ,t) = U,(Z,t) cos px + U,(Z,t) sin px, a9) 
v(x,Z,t) = V,(Z,t) cos ux + V,(Z,t) sin ux, (2.20) 
w(x,Z,t) = W,(Z,t) cos ux + W,(Z,t) sin px, ea, 
@(x,Z,t) = H,(Z,t) cos ux + H,(Z,t) sin wx, (e222) 
Q(x,Z,t) = Q,(Z,t) cos nx + Q,(Z,t) sin ux and Cpe) 
MTS(x,Z,t) = MTS,(Z,t) cos ux + MTS,(Z,t) sin ux. (2.24) 


The relations (2.15) through (2.24) are substituted into equations (2.5) through 
(2.12). The prognostic and diagnostic equations are separated into equations for the 


cosine and the sine terms. The resultant prognostic equations are 
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The resultant diagnostic equations for u and v are 


D, 


U,=- Tae (2733) 
5 BY 
L, = ose (2.34) 
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Geopotential values above the surface are obtained bv integrating the hydrostatic 


equation from the surface (Z = Z,) to height Z. 


Z 
H, = R| T,(Z,t)}dZ +S, and (2.37) 
Ze 


H,=R “TAZ, t)dZ + S,. (2.38) 
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The vertical velocity is calculated by integrating the continuity equation froin the top 
of the atmosphere (Z = Z,) down to height Z. The upper boundarv condition, w = 0 


at Z = Z,, is used. The diagnostic equations for the vertical velocity are 


Z 
Ww, -| 'D, (Z,t)dZ and (2559) 
z 


Té 
Wy, = | "D,(Z,t)dZ. (2.40) 
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Equations (2.25) through (2.40) are the prognostic and diagnostic equations that 
govern all six numerical models. Using the given basic state and the one-wave spectral 
perturbation quantities, the governing equations reduce to functions of Z and t. The 
models are effectively one-dimensional (Shapiro, 1987). 

To display the results of each model, the sine and cosine amplitudes of each variable 
are combined to determine the amplitude and phase of a single cosine wave in the x- 


direction. A typical variable has the form: 
mix. Z t= A(Z, ti costux — 0), (2.41) 


where the amplitude is A(Z,t) and the phase is 6(Z,t). The amplitude and phase are 


calculated at each level for all variables. 


C. TIME DIFFERENCING 

Two forward time steps are taken to start each model and then leapfrog time dif- 
ferencing 1s used. The leapfrog scheme is a higher order scheme than forward differ- 
encing and is the one used in many numerical models. A Robert filter is used to reduce 
the amplitude of the computational mode generated by the leapfrog time differencing. 
The filter is discussed by Haltiner and Wilhams (1980). For a prognostic variable F, 


calculate Fbar,_,, the average value of F at time step (n-1)At, using equation (2.42), 
papere tira er. |e re) (2.42) 


where y is a weighting function. Using the unaveraged values at time step nAt, compute 
the tendency (¢F/ét), from its predictive equation. The predicted value at time step 


(n+ 1)At is then calculated using equation (2.43), 


Fyat = Fao + 24t( 5), (2.43) 


In all the experiments, y = 0.05 is used. The time step for each experiment is calculated 


in the model by requiring, for computational stability, 


vAt = 


rl 


: (2.44) 


where v = uc and Cc is the typical phase speed of an external gravity wave. 


D. VERTICAL GRIDS 

Each of the models uses one of three vertical grids. The three methods for distrib- 
uting the variables over discrete levels are depicted in Fig. 1. The staggered levels are 
represented by the dashed lines in Fig. 1. Notice that the heights at which the variables 
are defined change between the three grids. The notation used in this paper to denote 
the staggered and unstaggered levels is consistent with the conventions used in the coded 
models. The height of the unstaggered levels 1s denoted as Z’. The height of the stag- 
gered levels is denoted as Z. In the models, both Z,’ and Z, are defined to be the surface 
of the earth. It is assumed that the staggered level Z; is exactly in the middle of the laver 
between Z;_,’ and Z,;’. This distinction 1s important because the models can have layers 
with unequal depth. Thus, the height of the staggered levels is defined relative to the 
height of the unstaggered Icvels. 

A finite difference model is written for each of the grid structures. The models are 
denoted as FDM-A, FDM-B and FDM-C. Similarly, finite element models using the 
three grids are indicated by FEM-A, FEM-B and FEM-C. 


E. FINITE DIFFERENCE MODELS 

The only differences in the equations between the three FDM models are the ap- 
proximations of terms involving dubar/dZ and éTbar/dZ in the prognostic equations and 
the approximations of the integral in the diagnostic geopotential equation. Centered 
difference approximations are used, except at the boundaries where one-sided differences 
are used. The finite difference approximations used in the prognostic equations are 


listed in Appendix A. 


F. FINITE ELEMENT MODELS 
]. FEM-A 
The FEM-A model defines vertical velocity (w) at the unstaggered levels in 
terms of the basis functions w,(Z). The other variables are defined at the staggered levels 


in terms of the basis functions @(Z). The expansion for a typical term is 


n+1 n+1 


AZ) = ) AMOW(Z) or ) Aloo) (2.45) 


The basis functions for this model are depicted in Fig. 2. The basis functions w(Z) are 
defined for the unstaggered levels (solid lines at height Z’) and the basis functions @(Z) 


are defined for the staggered levels (dashed lines at height Z). 





Fig. 2. Basis functions for grids A and B. 
The finite element approximations for the vorticity, divergence and 
thermodynamic equations are derived by substituting the expansion for each dependent 
variable into equations (2.25) through (2.30). Each equation is multiplied by ¥(Z) and 


integrated with respect to Z from the bottom to the top of the atmosphere. Each term 
in the equations is the finite sum of separate integrals. Only the integrals of overlapping 
basis functions are nonzero. The resultant equations, listed in Appendix B, are matrix 
equations. For an n-layer model, the vectors A,, A,, D,, D., H,, H2, Q), QO) fe tee 
U,. V, and V,, contain n+2 components. The vectors W, and W, contain n+1 com- 
ponents. The matrices M, K, and ®, defined below, are (n + 2) x (n+ 2) matrices. The 
matrix P, defined below, is an (n+1)x(n+2) matrix. The mass matrix M for this 


model is defined by 


Zs 
M,, = $ (Z)b(Z)dZ for li-jl <1. (2.46) 


The matrix K is defined for terms multiplied by u, 


K;(u) = ya ‘| es (Z)O4(Z)b(Z)dZ for ji-—j| <1. (2.47) 
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The matrix P is defined for terms multiplied by ss W, or by a 


ar G 
P(X) = y x* cas Ws; (Z)b (Z)dZ for ji—j| < J, (2.48) 
kai~} 40 


where x 1s u or [. 


The matrix ® is defined for terms multiplied by ora N, 
aa 
®,(u) = > “|e $(Z)b(Z)dZ for ji—jl <1. (2.49) 


The staggered basis functions present two general problems for evaluating the 
elements of the four matrices. First, for an n-layer model, portions of basis functions 
p,(Z) and @,_,(Z) are defined in the model atmosphere but the physical meaning of 


contributions from those terms is unclear. The contributions are included in the first 


two rows and the last row of each matrix. Second, only portions of basis functions 
~{Z) and @,.,(Z) are defined in the model atmosphere. To describe the incomplete sides 
of both basis functions an assumption must be made about the value of ¢, at the surface 
and @,., at the top of the atmosphere. 

Assumptions are made and procedures are developed in an attempt to resolve 
these problems. In this model, the mean state variables, ubar and Tbar, are defined only 
at the n staggered levels. However, ubar and Tbar values defined at the nodal points 
of $,(Z) and @,.,(Z) are important in the Galerkin formulation of the dubar/dZ and 
OTbar/oZ terms. In these experiments, the values of ubar and Tbar are defined at the 
surface and top of the atmosphere. Jordan (1985) did not define them at the nodal 
points of @,(Z) and ¢,.,(Z). One of Shapiro’s (1987) major modifications of these ex- 
periments was to define ubar and Tbar at the nodal points of ¢,(Z) and @¢,.,{Z). For 
constant shear with height, ubar and Tbar are defined at the boundaries such that the 
shear in the two half lavers at the boundaries is the same as the shear in the other lavers. 
To evaluate the staggered basis functions defined in the layers between the surface and 
Z,, and Z,., and the top of the atmosphere, it is assumed that the value at the boundaries 
of those basis functions is one-half. Thus, three-fourths of the basis functions 
@({Z) and @,.,(Z) are defined in the model atmosphere. 

The equations for the general elements of the four matrices are evaluated by 
substituting into equations (2.46) through (2.49) the formulas for ¢,.,(Z), 6(Z). @,_,(Z), 
fen), w(Z), w.(Z), and y,_(Z), mm terms of the local coordinate €=Z—Z. The 
equations for these basis functions defined for levels 1, 2, 1, and n+1 are listed in Ap- 
pendix C. The matrices were evaluated by integrating numerically using 2 point 
Gaussian Quadratures as in Shapiro (1987). 

The vorticity, divergence and thermodynamic equations, written in matrix and 


vector form, are 





dA 

M —— = M(~fD, — BV,) — 2K(B)Az, (2.50) 
dA, 

M —— = M(-fD, — BV2) + wK()Ay, (2.51) 
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dD, 

M —+ = M(fA, — U2 + w°H,) + uK(G)D, — uP(@)W,, (2.53) 
dT, We Loa a 

M —— = — uK(O)T, + > SMV, — PTW, + MQ, and (2.54) 
dT, ciel f = 

M te uK(u)T, + R @(u)V, — P(T)W, + MQ. (2755) 


Equations (2.50) through (2.55) are simplified by multiplying each equation by 
M-' and applying the Robert filter. Actually, one should not compute the inverse of 
M. Instead, at t=O one should obtain the LU factorization of M. Thus, at each time 
step one onlv needs to forward and back solve a triangular system. The matrices 
M-'K, M-'P and M-'® are constants. They are constructed in the initialization sub- 
routine and stored for use in the forecast subroutine. The matrices are multiplied by the 
appropriate vectors with values for time level nAt. The resultant forecast equations are 
vector equations and the forecast value for the i-th vertical level is the sum of values in 
the i-th location of each vector equation. The prognostic equations for the vorticity, 


divergence and potential temperature vectors are 

Aynsty = Arner + 2At( -fD, — BV, — EM 'K()Ag)iny (2.56) 
Anni) = Agn-y + 2At( —fD, — BV2 + EM™'K(Q)A) ony, (2.57) 
Dine = Di (m1) 

+ 2A(fA, — BU, + 2H, — »M7'K(@)D, — uM™'P(Q)W,)q, (Qs) 
Dan41) = Da(n=1) 

+ 2AU(fA, — BU + pH, — uM7'K(U)D, + uM7'P)W) ay, (2.59) 
Tatty = Tim-1 


+ 2At(— uM~'K()T, + + M7 'O@)V, — M“'P(T)W, + Qy)jq) and (2.60) 
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Tangy ie Tan ) 
+ 2At(uMK(B)T, + + M7 '@(G)V, — M7 P(T)W, + Q5)ep (2.61) 


where the subscripts (n+ 1), (n) and (n-1) refer to the values of the vectors at time step 
(n+1)At, nAt and (n-1)At, respectively. The surface geopotential and the diagnostic 
variables are calculated using the corresponding equations in model FDM.-A (see Ap- 
pendix A). 

2. FEM-B 

The FEM-B model defines vertical velocity, potential temperature, mean state 
potential temperature and diabatic heating at the unstaggered levels in terms of the basis 
functions y(Z). The other variables are defined at the staggered levels in terms of the 
basis functions @(Z). The basis functions are the same as defined for the FEM-A model, 
shown in Fig. 2. 

The finite element approximations for the vorticity, divergence and 
thermodynamic equations are derived by substituting the expansion for each dependent 
variable into equations (2.25) through (2.30). The vorticity and divergence equations 
are multiplied by @(Z) and integrated with respect to Z from the bottom to the top of 
the atmosphere. The resultant Galerkin formulation of the vorticity and divergence 
equations are the same as those derived for model FEM-A. The matrices in those 
equations, M, K and P are given by equations (2.46) through (2.48) as defined for 
FEM-A. The thermodynamic equations are multiplied by wZ) because potential tem- 
perature is defined at the unstaggered levels. As before, the equations are integrated 
through the depth of the atmosphere. The resultant equations are listed in Appendix 
D. Four additional matrices are defined for the two thermodynamic equations. The 


mass matrix II is 


zs | - 
Ml, = I (ZW ZIAZ for |i—jl <1. (2.62) 


The matrix Tis defined for terms multiplied by u , 


i+] 





Ze an 
P(u) = : a| W(Z)b,(Z)w(Z)dZ for |1—j| < 1. (2.6359 
oe 
Terms multiplied bv ria V, give rise to the transpose of the matrix P, defined by (2.48). 
The matrix W is defined for terms multiplied by ora W, 
i+] : 
= —,{ | dy, oe 
k=1—1 LZ 


As discussed in the FEM-A model description, the staggered finite elements 
present problems for evaluating the elements of the matrices. In this model, ubar 1s de- 
fined at the surface, the top of the atmosphere and at the n staggered levels. The mean 
state temperature, Tbar, is defined at the unstaggered levels so special definitions for it 
are not needed. Jordan (1985) did not include the contributions from the perturbation 
quantities defined at the nodal points of ,(Z) and @¢,.,(Z). They were included in this 
model as part of the Shapiro (1987) modifications. The staggered basis functions, 
@(Z), are evaluated at the boundaries using the assumptions discussed in the previous 
section. 

The elements of matrices I, , and W are evaluated by substituting formulas 
for ¢,,..(Z), $,.,(Z), 6(Z), $,-,(Z), WZ), w{Z), and w_,(Z) defined in terms of the local 
coordinate € = Z - Z, into equations (2.62), (2.63) and (2.64). Formulas for these basis 
functions are listed in Appendix C. As in FEM-A, the matrices are evaluated by inte- 
grating numerically using 2 point Gaussian Quadratures. 

The forecast matrix equations for vorticity, divergence and temperature are as 
described for model FEM-A. The final form of the vorticity and divergence vector 
equations are the same as for model FEM-A, equations (2.56) through (2.59). The 


thermodynamic vector equations are 


Th (n41) = yess 


+ 2At( —uHT TT, + + NPV, — NN YD)W, + Qn, (2.65) 


Tongs) a Tyg ) 


+ 2At(ulI F(T, + > MP 'a)V, — 1H (T)W, + Qn) (2.66) 
Pn tiniisemodem tae vectors Ay, A, DD Hy, H,, U,, U,, V,, V., W, and W,, 
contain n+2 components. The vectors Q,, Q,, T, and T, contain n+1 components. 
Shapiro (1987) defined vectors W, and W, with n+ 1 components; however, since W, and 
W, are derived using vectors D, and D, (which have n+2 components) in equations 
(2.39) and (2.40), W, and W, have been increased to n+2 components. The matrices 
II, [, and WY, are (n+ 1)x(n+1) matrices and the matrix P? is an (n+ 2) x (n+ 1) 
matrix. The surface geopotential and the diagnostic variables are calculated using the 
corresponding equations in model FDM-B. 
3. FEM-C 
The unstaggered FEM model is the simplest of the three FEM models. Each 
of the dependent variables is expanded in terms of the basis functions @(Z). The ex- 


pansion for a typical term is 


n+] 
A,(Z,t) = ) Aloo) (2.67) 


j=! 


The basis functions for this model are depicted in Fig. 3. 

The finite element approximations for the vorticity, divergence and 
thermodynamic equations are derived by substituting the expansion for each dependent 
variable into equations (2.25) through (2.30). Each equation 1s multiplied by @(Z) and 
integrated with respect to Z from the bottom to the top of the atmosphere. Each term 
in the equations is the finite sum of separate integrals. Only the integrals of overlapping 
basis functions are nonzero. The resultant equations, listed in Appendix E, are matrix 
equations. The matrices M and K are given by equations (2.46) and (2.47) as defined for 


FEM-A. The matrix E 1s defined for terms multiphed by a V, by Ww, or by 
for 
a7 aE 





Fig. 3. Basis functions for grid C. 


i+1 
E _x{ "dd, a 
Ei(x) = x Tv &(Z)0(Z)dZ fOr [t= 4[-= 2 (2.68) 
k=i-)  * 40 


where X 18 u or T. 


Note that this is the same as equation (2.49) which 1s only defined for x =u and is an 
(n+ 2) x (n+ 2) matrix. It has been given a new name because it combines equations 
(2.48) and (2.49) and is now an (n+ 1) x (n+ 1) matrix. The matrices M, K and E are 
evaluated using the same procedures as models FEM-A and FEM-B except only the 
formulas for $,(¢) are needed from Appendix C. 

The vorticity and divergence equations, written in matrix and vector form, are 
given in equations (2.50) through (2.53) as defined for model FEM-A. The 


thermodynamic matrix equations are 


dT 7 
M = = — »K(a)T, + > E(u)V, — E(T)W, + MQ,, (2.69) 


dT, _ _ = 
M —> =uK@)T, + ~ E(u)V, — E(T)W, + MQ,. (2.70) 


These forecast matrix equations are simplified as described for model FEM-A. 
The final form of the vorticity and divergence vector equations are the same as for model 


FEM-A, equations (2.56) through (2.59). The thermodynamic vector equations are 


Tints) = eres 


+ 2At( — pM 'K(a)T, + + M7'E(a)V, — M7E(T)W, + Qy)q) and (2.71) 


Tangy = Ty-1) 
+ 2A(uM~'K(G)T, + M“E(@)V, — MTEC) W) + Qi.) (2.72) 


emeiceimodel tie vectors A,, A,, DD, Fly H., ©,, Q,, T,, I., by, Ga Ve V3, 
W, and W, contain n+l components. The matrices M, K and E are (n+ 1)x(n+1) 
matrices. The surface geopotential and the diagnostic variables are calculated using the 


corresponding equations in model FDM-C. 


Il. EXPERIMENTS AND RESULTS 


Several experiments are performed for baroclinic instability with all or some of the 
models. The first results use the linear wind profile which was used by Shapiro (1987). 
Most of the experiments use a hyperbolic tangent vertical wind profile with three vari- 
ations. The basic experiments use an initial v-field which is independent of Z with no 
temperature perturbation. These experiments are integrated for 96 hours to determine 
the exponential growth rate. The basic experiments use two-, four-, six- and sixty-layer 
models. The next to last experiment uses the finite difference method of calculating the 
thermodynamic variables in the finite element model for grid A. The last experiment 
uses an initial v-field which tilts backward (toward negative x) with height. The temper- 
ature 1s chosen to satisfy the thermal wind relation. These profiles are integrated out to 
72 hours to investigate the enhanced initial growth due to the tilt (Farrell, 1985). The 
last two experiments use only the sixty-layer models except for the thermodynamic ex- 
periment which includes six-laver models. The analytical solution of each experiment 
does not exist except for the Eady (1949) solution for the linear shear solutions. 

As in Jordan (1985) and Shapiro (1987), the sixty-laver model results are intercom- 
pared to determine if the models are converging to the same solution. The standard of 
comparison for the two-, four- and six-layer models is the consensus sixty-laver solution. 
The v-component and temperature profiles are examined in each experiment. Expo- 
nential growth rates are calculated to help compare the models. 

Several parameters are defined identically in each experiment. The vertical coordi- 
nate, Z, 1s defined between zero and one (1000-368mb) and the vertical levels are equally 
spaced. The Coriolis parameter is defined at 45 degrees latitude. The mean state po- 
tential temperature increases with height from its surface value of 310.0 degrees Kelvin. 
The vertical shear in the ubar field is defined in each experiment as a function of Z 
multiplied by STRGTH, where STRGTH defines the strength of the wind at the top of 
the atmosphere in m/s. In all experiments, STRGTH is defined as 40 m/s. There is no 
diabatic heat source and no mountain topography. Unless otherwise indicated, all 
forecast experiments are terminated at 96 hours. 
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A. LINEAR VERTICAL SHEAR IN THE UBAR FIELD 

The wind profile in these experiments is a linear function of Z, with ubar = 
(STRGTH)Z. Waves are generated in each model using an initial perturbation of v’ = 
5.0 m/s in the cosine term. Wavelengths (in x) of 4000, 3000 and 2000 km are used for 
each model for two, six and sixty levels. 

1. Sixty-Layer Models 

For an x-wavelength of 4000 km, the sixty-layer FEM-C and FDM-C models 
both converge to the same solution as the FDM-A and FDM-B models from Shapiro 
(1987) for the v-component amplitude and phase (Fig. 4) and the temperature amplitude 
(Fig. 5). All figures will be found at the end of the chapter. As in Shapiro (1987), the 
phase is defined between zero and 360 degrees and a discontinuity occurs in the phase 
profile if the phase passes through zero degrees. FDM-A, FDM-B, FDM-C and FEM-C 
represent the same physical solution, which 1s called the consensus solution. These sol- 
utions all have the same patterns which were derived by Eady (1949) for the quasi- 
geostrophic model. The v-field has a maximum amplitude on the boundaries and tilts 
backward with height. The FEM-B amplitudes in Fig. 4 and Fig. 5 are less than the 
consensus which indicates a smaller growth rate. The FEM-B temperature field (Fig. 
5) also has an oscillation near the surface similar to what was found by Shapiro even 
though an additional element is added to the vertical velocity. The FEM-A scheme has 
a slightly larger growth rate than the consensus. This suggests that the FEM-A and 
FEM-B models need further improvement to take full advantage of their possibilities. 

For an x-wavelength of 3000 km, the sixty-laver profiles for the v-component 
and temperature show the same results as the 4000 km case except that the magnitude 
of the amplitudes are smaller and the profiles appear to contain some oscillation (Fig. 
6). The reduced amplitude agrees with the Eady (1949) theory which gives smaller 
growth rates for this wavelength. 

At a smaller wavelength of 2000 km, the v-component and temperature ampli- 
tudes show very little growth at all but have several oscillations throughout the atmos- 
phere (Figs. 7-8). The finite difference models and FEM-C converge to the same 
consensus solution. FEM-A is very close to the consensus profile, but has a slightly 
smaller amplitude in the middle lavers of the temperature profile and a slightly larger 
amplitude in the v-component profile. FEM-B has a smaller amplitude than the con- 
sensus profile for the v-component and oscillates between having a smaller amplitude to 
having a larger amplitude in the temperature profile. The phase profiles show that all 


the models except FEM-A and FEM-B converge to the same phase solution (Fig. 7). 


FEM-A is very close to the consensus solution but changes from being slightiv greater 
in the lower layers to being slightly less in the upper lavers. FEM-B is smaller than the 
consensus in the lower layers and larger in the upper layers. The lack of growth in this 
case agrees with the Eady theory which gives only neutral modes for this scale. Farrell 
(1985) has shown that a given initial state can funnel energy into the neutral modes of 
the problem. In this case, the most likely neutral waves which appear in the solutions 
are the two quasi-geostrophic solutions and internal gravity waves with small vertical 
scales. 
2. Six-Layer Models 

For the 4000 km wavelength, all the models except FEM-B and FEM-A con- 
verge to the consensus for the v-component amplitude (Fig. 9) and all the models except 
FEM-B are converging to the consensus solution for the temperature amplitude (Fig. 
10). FEM-B has the same shape as the consensus profile but has a smaller amplitude 
for temperature and the v-component. FEM-A has a much smaller v-component than 
the consensus. The v-component phase profiles show that all the models except FEM-A 
and FEM-B converge to the consensus phase profile (Fig. 9). The six-laver FDM-B 
model is closest to the consensus. 

At 3000 km, the difference among the six-layer models and between the six-laver 
models and the consensus profile for v-component and temperature amplitudes becomes 
more pronounced (Figs. 11-12). In general, the finite difference models are closer than 
the finite element models to the consensus profile for schemes B and C. The opposite 
is true for grid A except in the lower lavers of the temperature amplitude profile. 
FEM-B has a much smaller amplitude than the consensus. The phase profiles show the 
models are still converging to the consensus phase profile except for FEM-B (Fig. 11). 

At 2000 km, the models appear to be a little chaotic (Figs. 13-14). However, 
all the models depict the overall shape of the consensus profile except for the oscillations 
in the temperature amplitude profiles. FDM-B appears to represent the consensus better 
than the other models although FEM-B gives the same surface temperature prediction 
as the consensus profile. In the v-component amplitude profiles, the finite element 
models are closer then the finite difference models to the consensus profile for grids A 
and C but not B. However, the surface prediction for FDM-C is closer than the pre- 
diction for FEM-C. FDM-B is closest to the consensus profile for the v-component 
amplitude. The phase profiles show that the models for scheme A are closer to the 
consensus phase profile than are the models for schemes B and C (Fig.13). Clearly, the 


small oscillations can not be reproduced in any six-layer scheme. 
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3. Two-Layer Models 

The curves for the two-layer models are not shown to save space. For an x- 
wavelength of 4000 km, FDM-C represents the shape of the consensus profile of tem- 
perature and v-component better than the other models, but it has the benefit of an 
extra prediction level for these variables. FEM-C gives the best surface prediction for 
both the temperature and v-component. In the v-component amplitude profiles, the fi- 
nite difference models are closer to the consensus amplitude than are the finite element 
models for schemes A and B. However, these two schemes do not have enough resol- 
ution to depict the shape of the consensus profile. The phase profiles for the v- 
component show general agreement among the models in representing the shape of the 
consensus. 

At 3000 km, FDM-C again represents the shape of the consensus amplitude 
profile for temperature and v-component well. FEM-A is also close to the magnitude 
of the consensus profile. The other models give poorer results. 

At 2000 km, the magnitude of FDM-C overshadows all the other models. 
FEM-C and FDM-A also show a tendency for growth. In the v-component phase 
profiles, all the A and B models have a negative tilt instead of the positive tilt of the 
consensus profile. FDM-C has a positive tilt in the lower atmosphere and FEM-C has 


a positive tilt in the upper atmosphere. 


B. NONLINEAR VERTICAL SHEAR IN THE UBAR FIELD 
The nonlinear vertical shear is created by using a hyperbolic tangent function of Z 
sinular to the one used by Satyamurty ef al (1982). In this experiment, ubar = 
(STRGTH) tanh(«Z+x*) where « and x are given different values to contro] the magni- 
tude and location of the shear. Fig. 15 illustrates the three ubar profiles used in this 
experiment which are referred to as Type I, Type II and Type III. Waves are generated 
in each model using an initial perturbation of v’ = 5.0 m/s in the cosine term. The x- 
wavelengths of 4000, 3000 and 2000 km are used for each model. 
1. Hyperbolic Shear Type I 
For this part of the experiment, « = 2.0 and x = 0.0. This increases the mag- 
nitude of the low level ubar shear and forces the largest gradient to occur near the sur- 
face. 
a. Sixty-Layer Models 
Table 1 contains the exponential growth rates for all the sixty-laver models. 


For an x-wavelength of 4000 km, the growth rates for Type I indicate that the sixty-layer 


a, 


finite difference models and FEM-C all converge to the same growth rate after 96 hours. 


This is similar to the results from the linear wind profile. 


Table 1 
Exponential Growth Rates ( x 10-5) 
Hyperbolic Shear 


60-layer Models 


po FD M-A| FDM-B| FDM-C| FEM-A | FEM-B | FEM-C| 
| | LC 


ipet 
4000 km 


3000 km 
2000 km 


4000 km 
3000 km 
2000 km 


ivoe Tt 


4000 km 
3000 km 
2000 km 





The growth rate for FEM-A is slightly larger and the growth rate for FEM-B 1s slightly 
smaller. FDM-A, FDM-B, FDM-C and FEM-C all converge to the same solution for 
the v-component and temperature amplitudes and represent the consensus solution. The 
FEM-A model has a slightly higher amplitude than the consensus profile and the shape 
is the same as the consensus profile. The FEM-B model has a lower temperature am- 
plitude than the consensus profile, but the shape is the same except for an oscillation in 
the lower layers similar to the linear shear case. The difference between FEM-B and the 
consensus profile 1s greater in the lower layers than in the upper layers. For the v- 
component amplitude, the amplitude for FEM-B is again lower than the consensus 
profile, but there is no oscillation in the lower layers of the atmosphere. The profiles for 
the v-component phase are the same for all six models. The profiles are not shown be- 


cause the consensus profile is included with the lower resolution curves. 
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For an x-wavelength of 3000 km, the growth rates show the same the results 
as 4000 km except the growth rates are greater. The v-component and temperature 
profiles also show the same results as 4000 km. Again, the amplitudes are larger for the 
smaller wavelength. The v-component phase profiles show that FDM-A, FDM-B, 
FDM-C and FEM-C are identical, but FEM-A and FEM-B have a slightly larger phase. 

For an x-wavelength of 2000 km, the growth rates show the same results 
as 4000 km except the growth rates are smaller. The v-component and temperature 
profiles also show the same results as 4000 km, but the amplitudes are smaller. The 
amplitude profiles also appear to contain some oscillation. The v-component phase 
profiles are identical for all the models except FEM-A which has a shghtly larger phase. 
This profile clearly has the instability shifted to smaller wavelengths than the linear 
profile due to the vertical wind shear being larger near the surface. 

b. Six-layer Models 

The profiles for the six- and two-layer models are not shown to save space. 
Table 2 contains the exponential growth rates for six-, four- and two-laver models for 
Type I. For an x-wavelength of 4000 km, the growth rates show that all the models 
except FEM-B are converging to the growth rate of the consensus profile. The v- 
component amplitude profiles support the convergence of the models to the consensus 
profile. The finite difference models are closer to the consensus solution than are the 
finite element models for grids A and B, but FEM-C is closer than FDM-C. All models 
represent the low level shear well. Although the growth rates show that both FDM-A 
and FEM-C are closest to the consensus, the amplitude profiles show that FEM-C is 
closest to the consensus profile. The temperature amplitude also shows the convergence 
of the models to the consensus solution except for FEM-B. In this case, the surface 
prediction for both FEM-B and FDM-B is considerably less than the consensus. The 
v-component phase profiles show a moderate amount of negative tlt in the lower half 
of the atmosphere corresponding to the low level shear region. All the models have the 
Same phase profile shape, but scheme C 1s closer to the consensus solution than are 
schemes A and B. This difference is due to a difference in phase speed between the 


schemes. 


2 


Table 2 
Exponential Growth Rates ( x 10) 
Hyperbolic Shear Type | 


Models 


|| Consensus | FDM-A| FDM-B| FDM-C| FEM-A | FEM-B | FEM-C_ 


6 Lavers 


4000 km 1.45 4: 1.42 1.49 1.43 1.36 1.46 
3000 km log : 1.48 1.54 1.49 1.41 1.38 
wins km a ce a ae a lead 


4 | 4 Layers | 
4000 km 
3000 km 
2000 km 


2 easers 


4000 km 
3000 km 
2000 km 
* These cases do not converge to a steady aes rate 





The growth rates for an x-wavelength of 3000 km give results similar to 
those for 4000 km except that FDM-C is closer than FEM-C to the consensus solution. 
The v-component amplitude and phase profiles are also similar. In this case, however, 
FDM-A 1s closest to the consensus amplitude profile instead of FEM-C, but FEM-B 1s 
still farthest from the consensus. The phase results are identical to the results for the 
4000 km wavelength case. 

At an x-wavelength of 2000 km, a few more changes are observed from the 
longer wavelengths. For this wavelength, the growth rates indicate that all the six-layer 
models are growing faster than the consensus solution. Although the growth rates for 
FDM-A and FEM-A are the same, the difference between the finite difference and finite 
element models for grids B and C is much greater at the smaller wavelength. These dif- 
ferences are also observed in the v-component and temperature amplitude profiles. Al- 
though the growth rates for FDM-A and FEM-A are the same, the amplitude profiles 


show that FEM-A 1s shghtly closer to the v-component consensus solution, but FDM-A 
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is slightly closer to the temperature amplitude consensus. This may occur because the 
growth rate is computed at the first laver above Z = O and the differences between 
FDM-A and FEM-A occur above this level. 

c. Four-Layer Models 

In the 4000 km wavelength case, the growth rates show that the four-layer 
models do not converge to the consensus solution as well as the six-layer 4000 km 
models. Otherwise, the results are basically the same. The finite difference models are 
closer than the finite element models for grids A and B, but not grid C. However, the 
v-component amplitude profiles indicate that FEM-C 1s closest only at the second layer, 
where the growth rate is computed, and the top laver (Fig. 16). The amplitude profiles 
also show the lesser convergence of the four-layer models to the consensus solution. 
The v-component phase profiles give the same results as the six-laver models, but with 
a slightly greater difference between the four-layer models and the consensus. 

For the 3000 km wavelength, the growth rates for the four-laver models 
again indicate results similar to the six-layer models. The finite difference models are 
closer to the consensus solution than are the finite element models, but FDM-B 1s clos- 
est to the consensus growth rate instead of FDM-A. FEM-A has the worst growth rate 
in this case. Also, the four-layer models do not converge to the consensus solution as 
well as do the six-laver models (Fig. 17). The v-component amplitude profiles support 
the growth rate results. The phase results for the four-laver models are similar to the 
results for the six-layer models. The primary difference is that except for FEM-C the 
four-laver model phase profiles are not as close to the consensus phase profile as are the 
six-layer models. 

For the 2000 km wavelength, the growth rates indicate much less conver- 
gence of the four-laver models to the consensus solution. In fact, FDM-C fails to con- 
verge to a steady growth rate altogether. For this case, the finite element model growth 
rates are closer than the finite difference model rates to the consensus solution for 
schemes A and C, but not scheme B. All the models except FEM-A and FDM-C grow 
much faster than the consensus. FEM-C is closest to the consensus growth rate. 
FEM-C also reflects the shape of the v-component amplitude consensus profile better 
than the other four-layer models (Fig. 18). All the models except FDM-C depict the 
same shape as the phase consensus and FEM-C is again closest to the consensus sol- 


ution. Note that the shorter waves also have smaller vertical scales. 
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d. Two-Layer Models 

For the 4000 km wavelength, the growth rates indicate that the two-laver 
finite difference models are closer to the consensus solution than are the finite element 
models for grids A and C, but not grid B. FEM-C 1s closest to the consensus amplitude 
profile for v-component and temperature. Schemes A and B do not have enough resol- 
ution to depict the low level shear in the v-component amplitude profile. The grid C 
models attempt to reflect the low level shear, but the surface predictions are considerably 
less than the consensus solution. All six models reflect the negative tilt of the consensus 
phase profile even though they have only two levels. FEM-C 1s closest to the consensus 
phase profile. 

The growth rates for the 3000 km case indicate that the finite difference 
models are closer than the finite element models to the consensus growth rate for grids 
A and B, but FEM-C 1s closer than FDM-C. The scheme B models and FEM-A have 
much slower growth rates than the consensus and FDM-C has a larger growth rate. 
FEM-C is closest to the consensus amplitude profile for the v-component, but FDM-A 
is closest in the middle and lower levels of the temperature profile. Although grid B also 
predicts the temperature at the surface, its prediction 1s much poorer than the consensus 
solution. The two-laver A and B models do not have enough resolution to depict the 
shape of consensus v-component phase profile but they still reflect the negative tilt. 

For the 2000 km wavelength, the growth rates for FDOM-A, FDM-C and 
FEM-C are much greater than the consensus growth rate, and FEM-A 1s slightly less 
than and closest to the consensus. FDM-B and FEM-B grow very little and their growth 
rates do not converge to a steady growth rate. The finite element models are closer than 
the finite difference models to the consensus v-component and temperature amplitude 
profiles for schemes A and C, but not for scheme B. FEM-A is closest in magnitude to 
the consensus solution. 

2. Hyperbolic Shear Type II 
For this part of the experiment, « = 4.0 and x = 0.0. This increases the mag- 
nitude of the low level ubar shear but keeps the largest gradient near the surface. 
a. Sixty-Layer Models 

For an x-wavelength of 4000 km, the growth rates for Type II indicate that 
the sixty-layer A and C models all converge to the same growth rate after 96 hours 
(Table 1). The growth rate for FDM-B is nearly the same but slighly smaller and the 
growth rate for FEM-B is the smallest. FDM-A, FDM-C, FEM-A and FEM-C all 


converge to the same solution for the v-component and temperature amplitudes and 
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represent the consensus solution. The FDM-B model has a slightly smaller amplitude 
than the consensus profile, but the shape is the same as the consensus profile. The 
FEM-B model has a lower temperature amplitude than the consensus profile and the 
Shape 1s the same except for a similar oscillation to that of the linear shear case in the 
lower layers. The difference between FEM-B and the consensus profile is greater in the 
lower layers than in the upper layers. For the v-component amplitude, the amplitude 
for FEM-B is again lower than the consensus profile, but there is no oscillation in the 
lower layers of the atmosphere. The magnitudes of the growth rates and the amplitude 
profiles are greater than for the Type I case. The profiles for the v-component phase 
are the same for all six models. The sixty-layer profiles are not shown because the con- 
sensus profile is included with the lower resolution curves. 

For an x-wavelength of 3000 km, the growth rates for all the sixty-layer 
models are identical except for FEM-B and are greater than the 4000 km growth rates. 
The v-component and temperature profiles show results similar to the 4000 km case, but 
with larger amplitudes. The v-component phase profiles show that all the models have 
identical phases at 96 hours except for FEM-A which has a slightly larger phase. 

For an x-wavelength of 2000 km, the results are the same as for the 3000 
km case except FEM-C has a slightly larger growth rate. However, the v-component 
and temperature profiles show the same results as for 4000 km except the amplitudes are 
larger for 2000 km. The v-component phase profiles are identical for all the models ex- 
cept FEM-A which has a slightly larger phase. 

b. Six-Layer Models 

To save space, the profiles for the six- and two-layer models are not shown. 
Table 3 contains the exponential growth rates for six-, four- and two-layer models for 
Type IJ. For an x-wavelength of 4000 km, the growth rates indicate that the finite ele- 
ment models are closer to the consensus growth rate than are the finite difference models 
for grids B and C, but not for grid A. The growth rates for grid C are closer than the 
rates for grids A and B to the consensus. FEM-C is closest to the consensus growth 
rate. All the models represent the shape of the consensus profile well except for the grid 
B temperature prediction. FEM-C is closest to the consensus amplitude profile. The 
v-component phase profiles show the negative tilt in the lower half of the atmosphere 
corresponding to the region of low level shear. All the six-layer models have the same 
shape as the consensus phase profile, but scheme C is closer to the consensus than are 


schemes A and B. 
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Table 3 
Exponential Growth Rates ( x 10°*) 
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For the 3000 km wavelength case, the growth rates indicate results similar 
to the 4000 km case. For grids A and C, the finite difference models are closer to the 
consensus profile than are the finite element models although the differences are greater 
in the low level shear region than in the upper levels. FEM-B is closer than FDM-B to 
the consensus amplitude profile except at the midpoint of the temperature profile. 
FDM-C 1s closest to the consensus profile except in the lowest levels of the temperature 
profile above the surface where FEM-B is closest. The phase profiles show the same 
results as the 4000 km case. For the 2000 km wavelength case, all the results are the 
Same as the results for 4000 km. 

c. Four-Layer Models 

The growth rates for the four-layer 4000 km wavelength models indicate the 

finite difference models are closer than the finite element models to the consensus growth 


rate. Schemes A and B have much smaller growth rates than both scheme C and the 
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consensus solution. FDM-C has the same growth rate as the consensus. Overall, the 
v-component and temperature amplitude profiles show that the finite difference models 
are closer to the consensus profile than are the finite element models (Figs. 19-20). 
However, for scheme B the finite element model is closer to the consensus amplitude 
Paomle for the temperature. For scheme C, the finite element model is closer to the 
consensus at the surface and at Z = 0.25 for the v-component profile. Otherwise, 
FDM-C is closest to the consensus profile. The v-component phase profiles show the 
negative tilt in the lower half of the atmosphere corresponding to the region of low level 
shear (Fig. 19). All the six-layer models have the same shape as the consensus phase 
profile, but scheme C is closer to the consensus than are schemes A and B. 

For the 3000 km wavelength, the growth rates for the four-layer models 1n- 
dicate results similar to the 4000 km case, except that FEM-C 1s closer than FDM-C to 
the consensus solution. FEM-C is also closest to the growth rate for the consensus. 
Overall, the finite difference models are closer than the finite element models to the 
consensus amplitude profile for schemes A and B, but not for scheme C in the v- 
component profiles (Fig. 21). The finite element models are closer than the finite dif- 
ference models for grids B and C, but not grid A for the temperature profile (Fig. 22). 
FEM-C is closest to the consensus except in the middle lavers of the temperature profile 
where FDM-A is closest. For the phase profiles, the results for this case are the same 
as the results for the 4000 and 3000 km cases (Fig. 21). 

Growth rates for the four-layer models in the 2000 km case indicate that the 
finite difference models are closer than the finite element models to the consensus sol- 
ution for grids A and C. FDM-B and FEM-B are nearly identical. FIDM-C its closest 
to the growth rate for the consensus solution. In this case, the finite element models are 
closer than the finite difference models to the consensus amplitude profile for scheme B 
in the lower half of the atmosphere and for scheme C, but not for scheme A except in 
the upper levels of the v-component amplitude profile (Fig. 23). FEM-C 1s closest to the 
consensus profile. 

d. TIwo-Layer Models 

For the 4000 km wavelength, the growth rates indicate that the two-layer 
finite difference models are closer to the consensus solution than are the finite element 
models for grids B and C, but not for grid A. FEM-C 1s slightly closer than FDM-A to 
the consensus. FDM-C is closest to the growth rate for the consensus solution, but 
FEM-C is closest overall to the consensus amplitude profile for v-component and tem- 


perature. Schemes A and B do not have enough resolution to depict the low level shear 
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in the v-component amplitude profile. The grid C models attempt to reflect the low level 
shear, but the surface predictions are considerably less than the consensus solution. All 
six models reflect the negative tilt of the v-component even though they have only two 
levels and FEM-C is closest to the consensus phase profile. 

For the 3000 km case, the growth rates show that the finite difference 
models are closer than the finite element models to the concensus for schemes A and 
C, but not for scheme B. FEM-B appears to be closer than FDM-B to the consensus 
growth rate, but FEM-B does not converge to a steady growth rate. The growth rate 
for FDM-C is closest to the consensus. The amplitude profiles for scheme C are also 
closest to the consensus amplitude profile for the v-component and temperature, and 
schemes A and B have much smaller magnitudes than the consensus. The finite differ- 
ence models are closer than the finite element models to the consensus for schemes A 
and C except for FEM-C at the top laver. FEM-B is closer than FDM-B. For the v- 
component phase, all of the models depict the negative tilt of the consensus profile. 

At 2000 km, the growth rates for grid C are greater than the consensus, and 
the growth rates for grids A and B are much less than the consensus growth rate. The 
growth rates for grid B do not converge to a steady growth rate. The growth rate for 
FEM-C is closest to the consensus. The finite element models are closer to the consensus 
amplitude profile for temperature and v-component than are the finite difference models 
for grids B and C, but not for grid A. Scheme C has a greater magnitude than the con- 
sensus, especially FDM-C. FDM-A 1s closest to the consensus in magnitude but does 
not have enough resolution to depict the low level shear. 

3. Hyperbolic Shear Type III 
For this part of the experiment, a = 4.0 and x = 1.0. This causes the strongest 
shear to occur at about Z = 0.25 instead of at the surface, and provides a variety of 
responses that are less sensitive to the lower boundary conditions. 
a. Sixty-Layer Models 

At the 4000 km wavelength, all the growth rates are identical except for 
FEM-B which has a smaller growth rate (Table 1). All of the sixty-layer models except 
FDM-B and FEM-B converge to the same solution for the v-component and temper- 
ature amplitude profiles. FDM-B is identical to the consensus profiles except at the 
surface of the temperature profile. FEM-B shows the same results as the Type II case, 
but the temperature oscillations are smaller (Fig. 24). The magnitudes of the growth 
rates and the amplitude profiles are larger than those for Type I and Type II. The phase 
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profiles of all the models converge to the same solution except FEM-B which has a 
slightly larger phase. 

For 3000 km, FDM-B, FDM-C and FEM-C have identical growth rates. 
FDM-A and FEM-A also have identical but slightly smaller growth rates and FEM-B 
has the smallest growth rate. The v-component and temperature amplitude profiles 
show small differences between the models, especially near the strongest shear region at 
Z = 0.25 (Figs. 25-26). For the v-component amplitude profiles, the finite difference 
models converge to the same solution and represent the consensus solution. FEM-A 
has a slightly smaller amplitude and FEM-C has a slightly larger amplitude 1n the shear 
zone. In the temperature profiles, the differences are slightly larger for FEM-A and 
FEM-C. The results for FEM-B are the same as the previous cases, but again the os- 
cillations are smaller. For the v-component phase profiles, all the models converge to 
the same phase except FEM-A and FEM-B which are slightly larger (Fig. 25). 

At the 2000 km wavelength, the growth rates show a little more variation 
than in the other cases. FDM-B and FDM-C are identical and represent the consensus 
growth rate. Although verv close to the consensus, the growth rate for FDM-A 1s 
slightly smaller and FEM-C 1s a little larger. FEM-A has a slightly smaller growth rate 
than FEM-A and FEM-B has the smallest growth rate. The v-component and temper- 
ature amplitudes also show more variations among the models in the region of the 
stongest shear than in previous sixty-laver cases (Figs 27-28). FDM-B and FDM-C are 
identical and represent the consensus profile. FDM-A is very close to the consensus, 
but slightlv smaller in magnitude. FEM-A has a smaller magnitude than the consensus 
and FDM-A. FEM-C has a larger magnitude than the consensus. FEM-B has a much 
smaller magnitude than the consensus. The differences among the models are slightly 
greater for the temperature profiles than the v-component profiles. However, for the 
phase profiles, the finite difference models and FEM-C converge to the same phase. 
FEM-A 1s slightly greater and FEM-B is slightly smaller (Fig. 27). 

b. Six-Layer Models 

To save space, the profiles for the six- and two-layer models are not shown. 
Table 4 contains the exponential growth rates the six-, four- and two-layer models for 
Type III. For a wavelength of 4000 km, the growth rates for the finite element models 
are closer to the consensus than are the growth rates for the finite difference models for 
schemes A and C, but not B. FDM-B is closest to the consensus growth rate. Overall, 
the finite element models are closer than the finite difference models to the consensus 


amplitude profile for v-component and temperature for grids A and C. For grid B, the 
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finite difference model is closer to the consensus than the finite element model except 
at the lowest level above the surface. In general, grid B is closest to the consensus am- 
plitude profile. For the phase profiles, all of the models are close to the consensus pro- 


file, but grid B 1s closest. 


Table 4 
Exponential Growth Rates ( x 10-5) 
Hyperbolic Shear Type III 
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For the 3000 km case, the growth rates indicate that all the models are at 
least 8 percent less than the consensus growth rate. The finite element models are closer 
than the finite different models to the consensus and FEM-B is closest. Overall, the fi- 
nite element models are closer than the finite difference models to the consensus ampli- 
tude profile for v-component and temperature. The only exceptions are at the top layer 
of the v-component amplitude profiles and at the midpoint of the temperature profiles 
where FDM-B is closer. In general, FEM-B is closest to the consensus. For the phase 
profiles, the A and B profiles are closer than the C profiles to the consensus, and 


FDM-A is closest to consensus phase profile for the v-component. 
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For a wavelength of 2000 km, the growth rates for the finite element models 
are closer than the growth rates for the finite difference models to the consensus growth 
rate with FEM-A being closest. The growth rates for the finite difference models and 
FEM-C are much less than the growth rate for the consensus. The growth rates for the 
C models are still converging at 96 hours, but even after 192 hours the growth rates in- 
crease only shghtly. The amplitude profiles for the v-component and temperature also 
show that the finite element models are closer than the finite difference models to the 
consensus. FEM-A is closest to the consensus amplitude profile. For the phase profiles, 
the scheme A models are closer to the consensus phase profile than are the scheme B 
and C models. 

c. Four-Layer Models 

At an x-wavelength of 4000 km, the growth rates for FDM-B and FEM-B 
are the same and closest to the consensus growth rate. For schemes A and C the growth 
rates for the finite element models are closer to the consensus than are the growth rates 
for the finite difference models. The growth rates for schemes A and C are much less 
than the consensus. The finite element models are closer than the finite difference 
models to the consensus amplitude profile for grids A and C (Fig. 29). For grid B, 
FDM-B is closer than FEM-B and 1s the closest overall to the consensus profile. For 
the phase profiles, all the models follow the shape of the consensus profile well (Fig. 29). 
FEM-C is closest to the consensus phase profile. 

For 3000 km, the growth rates for the finite element models are closer than 
the finite difference models to the consensus for schemes A and C, but not B. FDM-B 
is closest to the consensus and schemes A and C are much less than the consensus. For 
grids A and C, the finite element models are closer than the finite difference models to 
the consensus amplitude profile for the v-component and temperature (Figs. 30-31). For 
grid B, FDM-B 1s closer than FEM-B to the consensus profile except in the upper lavers 
of the v-component profile. Otherwise, FDM-B is closest to the consensus profile for 
this case. For the phase profiles, FEM-C is closest to consensus phase profile (Fig. 30). 

For 2000 km, the growth rates show that none of the four-layer model 
growth rates are close to the consensus. The growth rates for the finite difference models 
are closer than the finite element models to the consensus growth rate. The growth rates 
for the grid B models are much greater than the consensus and the growth rates for the 
grid A and C models are much less. The growth rate for FDM-B is closest to the con- 
sensus growth rate. The amplitude profiles for the v-component and temperature show 


the same results as the growth rates (Figs. 32-33). The phase profiles show that all the 
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models except FDM-C have the same shape as the consensus profile and that scheme 
C is closest to the consensus (Fig. 32). 
d. Two-Layer Models 

For the 4000 km wavelength, the growth rates indicate that none of the 
models grow as rapidly as the consensus model for this case. For grids A and C, the 
growth rates for the finite difference models are closer than the rates for the finite ele- 
ment models to the consensus growth rate. FEM-B is closer to the consensus than 
FDM-B, and FDM-C is closest to the consensus. The finite element profiles are closer 
than the finite difference profiles to the consensus for scheme B and equal to or slightly 
closer to the consensus for scheme C. FDM-A is closer than FEM-A. FEM-Bi11s closest 
to the consensus in this case except at the surface for the temperature profile where 
FEM-C is better. However, all the surface predictions are considerably less than the 
consensus prediction. For the phase profiles, all six models reflect the negative tilt of the 
v-component even though they have only two levels. FEM-C 1s closest to the consensus 
phase profile. 

For the 3000 km case, the growth rates indicate the same results as the 4000 
km case. The amplitude profiles show that except for grid B, the finite difference models 
are closer than the finite element models to the v-component and temperature amplitude 
consensus profile. FDM-C is closest to the consensus profile. For the phase profiles, 
all the models follow the shape of the consensus phase profile except FDM-C. FEM-A 
is closest to the consensus phase profile for this case. 

At the 2000 km x-wavelength, the growth rates indicate that the grid C 
models have a much greater growth rate than the consensus model and that FEM-A and 
the grid B models have a much smaller growth rate. FDM-A also has a smaller growth 
rate than the consensus, but it is closer to the consensus than are the other models. 
Overall, the finite element models are closer than the finite difference models to the 
consensus amplitude profile for grids B and C, but not grid A. FDM-A is closest to the 
consensus profile in magnitude, but does not have enough resolution to depict the low 
level shear. For the v-component phase profiles, the models have a negative tilt similar 
to the consensus profile, but do not have enough resolution to depict the shape of the 


consensus profile. FDM-C is closest to the consensus phase profile. 


C. USING FINITE DIFFERENCE METHOD IN FINITE ELEMENT MODEL 
In this experiment, the finite difference equations and coefficients for the 


thermodynamic predictive equations are used in the finite element model for scheme A 
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instead of the matrix equations. The model is run with linear vertical shear in the ubar 
field for sixty-layer models and with hyperbolic Type III vertical shear in the ubar field 
for six- and sixty-laver models. Waves are generated in the model using an initial per- 
turbation of v = 5.0 m/s in the cosine term. Horizontal wavelengths of 4000 and 2000 
km are used. The results are compared with the unmodified FDOM-A and FEM-A pro- 
files. 
1. Linear Vertical Shear in the ubar Field 
For both 4000 and 2000 km wavelengths, the temperature and v-component 
amplitude profile of the modified FEM-A model 1s identical to the unmodified FEM-A 
model (Figs.34-37). The phase profiles show the same results. 
2. Hyperbolic Type III Shear in the ubar Field 
a. Sixty-Layer Model 
When the x-wavelength is 4000 km, the v-component amplitude shows that 
the modified sixtv-layer FEM-C model has a greater amplitude than the unmodified 
models throughout the entire atmosphere (Fig. 38). The temperature amplitude of the 
modified FEM-C model is greater than the amplitudes of the unmodified models in the 
region of greatest vertical shear in the ubar field (Fig. 39). The phase profiles show that 
the phase for the modified model is identical to the unmodified FEM-A model (Fig. 
38-39). For the 2000 km wavelength, the results are the same, but the difference between 
the amplitudes for the modified and unmodified models is greater (Fig. 40). 
b. Six-Layer Model 
For the 4000 km wavelength, the temperature and v-component amplitudes 
of the modified FEM-A model are nearly the same as the profiles for FDM-A (Figs. 
41-42). The phase profiles show the modified FEM-C model is closer to the unmodified 
FEM-C model except in the middle lavers of the temperature phase profile (Figs. 41-42). 
When the x-wavelength 1s 2000 km, the temperature and v-component am- 
plitude profiles show that the modified model lies between FDM-A and the unmodified 
FEM-C but is slightly closer to the finite difference model (Figs. 43-44). The phase 
profile of the modified FEM-C is closer to the unmodified FEM-C than to the finite 


difference model. 


D. BACKWARD TILTING INITIAL VALUES 
For this experiment, the initial phase in equation (2.41) is set to give a negative tilt by 
letting 6 = bZ for the vorticity and temperature initial conditions. The experiment is 


run for the sixty-layer FDM-C model with linear ubar shear and with hyperbolic Type 
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II] ubar shear. Each case is run with b = z and b = 4z which are called phase I and 
phase IJ respectively. The x-wavelengths of 4000, 3000 and 2000 km are used in each 
case. 
1. Linear Vertical Shear in the ubar Field 
a. X-wavelength of 4000 km 

The growth rates for this case do not converge to a steady growth rate al- 
though they appear to be converging for the 3000 and 4000 km wavelengths for phase 
I (Table 5). The phase profile for phase I shows a negative tilt of about 110 degrees at 
12 hours which becomes less negative by 72 hours (Fig. 45). For phase II, the phase 
profile has a much stronger negative tilt at I2 hours which also becomes more vertical 


with time. 


Table 5 
Exponential Growth Rates ( x 10-5) 
Linear Shear (Negative Initial Phase) 
Sixty-layer FDM-C Model 


X-wavelength (km) 
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b. X-wavelength of 3000 km 
The growth rates show the same general results as the 4000 km results ex- 
cept the growth rates are smaller for phase I and slightly greater until 36 hours for phase 
II. The phase profiles give the same results as the 4000 km case (Fig. 46). 
c. X-wavelength of 2000km 
In this case, the growth rates show no apparent convergence to a steady 
growth rate. The amplitude profiles show there 1s growth for about 24 hours and then 
the profile oscillates back and forth without growing (Fig. 47). This is expected because 
the Eady solution is not unstable for this wavelength. Phase II has the same pattern 
except it increases in magnitude for 48 hours before decreasing slighltv. For the phase 
profiles, phase I has a slightly negative tilt at 12 hours which becomes positive by 72 
hours (Fig. 48). Phase II has a greater negative tilt than phase | at 12 hours, but it is 
less than the tilt for the 4000 and 3000 km cases. The tilt for phase II becomes positive 
by 48 hours and then becomes negative again at 72 hours. 
2. Hyperbolic Type III Shear in the ubar Field 
a. X-wavelength of 4000 km 
The growth rates for the first 12 hours are much greater than the growth 
rates for the linear case, but converge quickly to a steadv growth rate for phase I (Table 
6). For phase II, the early growth rates are larger and converge slower than the rates 
for phase I. However, the growth rates for both phase I and phase IJ converge to the 
same value. The phase profiles show less negative tilt than the linear case (Fig. 49). 
Phase II has a larger negative tilt than phase I at 12 hours in the region of the strongest 
vertical shear in the ubar field (Z = 0.25). This negative tilt becomes less negative with 
time but the change is less pronounced than in the linear case. 
b. X-wavelength of 3000 km 
In this case, the results are similar to the results for 4000 km. The growth 
rates are greater than at 4000 km. The phase results are also the same as the 4000 km 


case except for less negative tilt at 12 hours (Fig. 50). 
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Table 6 
Exponential Growth Rates ( x 10°‘) 
Hyperbolic Shear Type III (Negative Initial Phase) 
Sixty-laver FDM-C Model 
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c. <A-wavelength of 2000 km 

At 2000 km, the same trend is still evident although the growth rates de- 
crease slightly from the 3000 km case after 12 hours. The growth rates do not oscillate 
as they did in the linear case. The phase profiles show the negative tilt at 12 hours is less 
than the previous cases (Fig. 51). 

The results of this experiment agree with the findings of Farrell (1985) that 
large initial growth rates can occur in the early time steps and converge to the most 
unstable growth. In the future, this type of time evolution should be used to test the 


various schemes for various wavelengths. 
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Fig. 4. 60-layer, 4000 km linear shear experiment at 96 h comparing v-component 
amplitudes (top) and phases (bottom). 
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Fig. 5. 60-layer, 4000 km linear shear experiment at 96 h comparing temperature 
amplitudes. 
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Fig. 6. 60-layer, 3000 km linear shear experiment at 96 h comparing v-component 
(top) and temperature (bottom) amplitudes. 
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LINEAR SHEAR (60-L, 2000KM) 
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Fig. 7. 60-laver., 2000 km linear shear experiment at 96 h comparing v-component 
amplitudes (top) and phases (bottom). 
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Fig. 8. 60-laver, 2000 km linear shear experiment at 96 h comparing temperature 
amplitudes. 
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Fig. 9. 6-laver, 4000 km linear shear experiment at 96 h comparing v-component 
amplitudes (top) and phases (bottom). 
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Fig. 10.  6-laver, 4000 km linear shear experiment at 96 h comparing temperature 
amplitudes. 
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Fig. 11. 6-layer, 3000 km linear shear experiment at 96 h comparing v-component 
amplitudes (top) and phases (bottom). 
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Fig. 12. 6-laver., 3000 km linear shear experiment at 96 h comparing temperature 
amplitudes. 
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13. 6-layer, 2000 km linear shear experiment at 96 h comparing v-component 
amplitudes (top) and phases (botton)). 
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Fig. 14. 6-laver. 2000 km linear shear experiment at 96 h comparing temperature 
amplitudes. 
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Fig. 16. 4-laver, 4000 km Type I shear experiment at 96 h comparing v-component 
amplitudes (top) and phases (bottom). 
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Fig. 17. 4-laver. 3000 km Type I shear experiment at 96 h comparing v-component 
amplitudes (top) and phases (bottom). 
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Fig. 18. 4-layver. 2000 km Type I shear experiment at 96 h comparing v-component 
amplitudes (top) and phases (bottom). 
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Fig. 19. 4-laver, 4000 km Type II shear experiment at 96 h comparing v-component 
amplitudes (top) and phases (bottom). 
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Fig. 20.  4-laver. 4000 km Type II shear experiment at 96 h comparing temperature 
amplitudes. 
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Fig. 21.  4-laver, 3000 km Type II shear experiment at 96h ere v-component 
amplitudes (top) and phases (bottom). 
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Fig. 22.  4-laver, 3000 km Type II shear experiment at 96 h comparing temperature 
anrplitudes. 
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Fig. 23. 4-laver, 2000 km Type IJ shear experiment at 96 h comparing v-component 
amplitudes. 
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Fig. 24.  60-laver, 4000 km Type II] shear experiment at 96 h comparing temper- 
ature amplitudes. 
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Fig. 25. 60-laver, 3000 km Type III shear experiment at 96 h comparing v- 
component amplitudes (top) and phases (bottom). 
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Fig. 26.  60-laver. 3000 km Type III shear experiment at 96 h comparing temper- 
ature amplitudes. 
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Fig. 27.  60-layer, 2000 km Type III shear experiment at 96 h comparing v- 
component amplitudes (top) and phases (bottom). 
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Fig. 28. 60-layer, 2000 km Type III shear experiment at 96 h comparing temper- 
ature amplitudes. 
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29, 4-laver, 4000 km Type III shear experiment at 96 h comparing v- 
component amplitudes (top) and phases (bottom). 
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Fig. 30. 4-layer, 3000 km Type II] shear experiment at 96 h comparing vy- 
component amplitudes (top) and phases (bottom). 
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Fie. 31. 4-laver, 3000 km Type III shear experiment at 96 h comparing temper- 
ature amplitudes. 
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Fig. 33.  4-laver, 2000 km Type III shear experiment at 96 h comparing temper- 
ature amplitudes. 
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Fig. 34.  60-laver, 4000 km linear thermodynamic modification at 96 h comparing 
v-component amplitudes (top) and phases (bottom). 
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Fig. 35.  60-laver, 4000 km linear thermodynamic modification at 96 h comparing 
temperature amplitudes (top) and phases (bottom). 
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Fig. 36. 60-laver, 2000 km linear thermodynamic modification at 96 h comparing 
v-component amplitudes (top) and phases (bottom). 
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Fig. 37. 60-layer, 2000 km linear thermodynamic modification at 96 h comparing 
temperature amplitudes (top) and phases (bottom). 
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Fig. 38. 60-laver, 4000 km Type II] thermodynamic modification at 96 h comparing 
V-component amplitudes (top) and phases (bottom). 
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Fig. 39. 60-layer, 4000 km Type III thermodynamic modification at 96 h comparing 
temperature amplitudes (top) and phases (hottom). 
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Fig. 41. 6-laver, 4000 km Type II] thermodynamic modification at 96 h comparing 
¥-component amplitudes (top) and phases (bottom). 
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Fig. 42. 6-layver, 4000 km Type II] thermodynamic modification at 96 h comparing 
temperature amplitudes (top) and phases (bottom). 
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Fig. 43. 6-laver, 2000 km Type III thermodynamic modification at 96 h comparing 
v-component amplitudes (top) and phases (bottom). 
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Fig. 44. 6-laver, 2000 km Type III thermodynamic modification at 96 h comparing 
temperature amplitudes (top) and phases (bottom). 
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Fig. 45. 60-laver, 4000 km initial value modification (linear) comparing phase I (top) 
and phase II (bottom) negative phase tilt. 
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Fig. 46. 60-laver, 3000 km initial value modification (linear) comparing phase I (top) 
_ and phase II (bottom) negative phase tilt. 
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Fig. 47. 60-laver, 2000 km initial value modification (linear) comparing phase I (top) 
and phase 1] (bottom) amplitudes. 
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Fig. 48. 60-laver, 2000 km initial value modification (linear) comparing phase I (top) 
and phase II (bottom) negative phase tilt. 
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Fig. 49.  60-layer. 4000 km initial value modification (Type III) comparing phase I 
(top) and phase II (bottom) negative phase tilt. 
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Fig. 50. 60-layver, 3000 km initial value modification (Type III) comparing phase I 
(top) and phase II (bottom) negative phase tilt. 
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Fig. S51. 60-laver, 2000 km initial value modification (Type II1) comparing phase 1 
(top) and phase II (bottom) negative phase tilt. 


IV. CONCLUSIONS 


With a linear vertical ubar profile, FEM-B still has an oscillation near the lower 
boundary of the temperature profile after correcting the length of the vertical velocity 
vector and is the poorest of the six models. FEM-A is better than FEM-B but is still 
not equal to the consensus profile. In the hyperbolic tangent vertical ubar profiles, the 
finite element models do not perform as well as the finite difference models for Types I 
and II except for the six-layer model for Type II. Grid A is closest to the consensus for 
most of the Type I and Type II cases. However, for Type III, the finite element models 
are better than the finite difference models in many of the cases. FEM-B shows much 
improvement in Type III. This may be because the maximum vertical shear in the ubar 
field occurs near Z = 0.25, which may be far enough from the boundary to enable 
FEM-B to overcome its earlier problems near the boundary. 

Overall, scheme C is often better than A and B, especially FEM-C. Scheme C mav 
have an advantage because it predicts the variables at one more level than the A and B 
models except for the vertical velocity and temperature. However, even for the temper- 
ature profile, where grid B predicts the same number of levels as grid C, scheme C often 
has the better surface prediction. Grid A seems to perform surprisingly well in spite of 
the analysis of Arakawa (1984). 

Future research should concentrate on the boundary elements for FEM-A and 
FEM-B. Perhaps new hybrid elements could be employed. Another approach would 
be to continue the work of isolating the effects of the vorticity, divergence and 
thermodynamic terms in the finite element models. This could be accomplished by using 
the finite difference method for these equations in the finite element models and ob- 
serving each term’s contribution. 

So far, all the baroclinic instability experiments have used a Tbar field that increases 
linearly with height, which may not be too challenging for simple finite difference models 
to solve. A vertical Tbar field that is not linear, such as the hyperbolic tangent profile 
used for the ubar field in this study, may allow the finite element models to take full 
advantage of their possibilities. Finally, future work should continue the initial value 
growth experiments by evaluating all the models and bv using smaller resolution models 
and smaller wavelength in x. Smaller horizontal wavelengths may lead to a different tvpe 


of growth. 
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APPENDIX A. FINITE DIFFERENCE APPROXIMATIONS 
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APPENDIX B. GALERKIN FORM OF FEM-A PROGNOSTIC 
EQUATIONS 


1. Vorticity Equations (2.25) and (2.26): 
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3. Thermodynamic Equations (2.29) and (2.30): 
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APPENDIX C. BASIS FUNCTION EQUATIONS FOR FENI-A 
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APPENDIX D. GALERKIN FORM OF FEM-B PROGNOSTIC 


EQUATIONS 


1. The vorticity equations have the same form as the vorticity equations for model 


FEM-A, equations (B.1) and (B.2). 


2. The divergence equations are the same as the divergence equations in model FEM-A, 


equations (B.3) and (B.4). 


3. Thermodynamic Equations (2.29) and (2.30): 


j=1—1 





Mi 


dT} 
pat | waren) Saf VibabndZ 
k=i-1  j=i—1 


ya)” AZ. 
a 


(D.1) 


ar enema) FY if V bur aZ 


k=1-1  j=l—] 


f od doy 
tys yu" p\wdZ 


=i~1  j=i—1 


ay ry) [S cue WiWidZ + dal WW,dZ. (D.2) 
k=i-1 j= t—-1 


Note that in these equations the basis functions are functions of Z (¢, = ,(Z) and 
w,=w(Z)) . All of the other variables, Q, T, u, U, V, and W, are functions of time 
(Q, = Q(t), T, = T(t), u,= u(t), U; = Ut), V, = V(t), W, = W(t) . 
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APPENDIX E. GALERKIN FORM OF FEM-C PROGNOSTIC 
EQUATIONS 


1. The vorticity equations have the same form as the vorticity equations for model 
FEM-A, equations (B.1) and (B.2). 


2. Divergence Equations (2.27) and (2.28): 
i+] i+] 


a |, ae = f ale ?,p,dZ 
dt 


j=1-1 j=i-1 


7S ‘yp Di] died eZ - »y vif PPA 


k=1-1 
i+] 


_ ys =k “le >) Hi "9,602, (E.1) 
Te 


j=1-1 j=1-1 
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+1 i+] 


>= $,,dZ = af $,b,dZ 


j=1-1 }=1-1 


Dy ye Di] didbeZ - ox vi] 8907 


j=i-1 
i+ 


Zt g BG 
yr SY wif fe gears YH |, ebaz (E.2) 
k=1-1  j=i-] 


j=i-1 





Note that in these equations, and the equations that follow, the basis functions are 
functions of Z (6, = ¢(Z) and ¥, = »(Z)). All of the other variables, A, D, H, Q, T, u, 
U, V and W, are functions of time (A;= A(t), D,= Dt), H,=H,t), Q, = Q(@), 
T, = T,(0, u, = u(t), 0, = b(t), V, = Vt) and W, = W((t)) . 


3. Thermodynamic Equations (2.29) and (2.30): 





Y at! 1+1 
ei epaze—n 1 Dn b6,,4Z 
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1+] i+] 7 
; Tq 
+ yey vil OO 6,42 
k=1-1 j=i-1 LZ 
i+] i+] 


— ey) e Gye — $.6,dZ + Yai)” 0 ;aZ. (E.3) 


k=1—-1 j=1--1 j=1-1 
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